## Shannon's limit and two Shannon points

from PyM import *

log10=np.log10

# Main functions
def shannon_limit(s): return dB((2**s-1)/s)
def shannon_point(s): return (shannon_limit(s),s)
def shannon_deficit(r,s): return r-shannon_limit(s) 

# Auxiliary functions
def dB(X): return 10*log10(X)
def dBi(x): return 10**(x/10)

# Convinient steps for the axis
dx=dy=0.1

# Important points
def X(k): return (k,0)
def Y(k): return (0,k)
B = (9.6,1); R = (6.26,1/3)
L = (-1.59,0); M=(-1.59,6+5*dy)

# Range for spectral efficiency (vertical axis)
# starts at 0.001 to avoid divisin by 0
s = ls(0.001,6,150)   

# plotting
close('all')

ax = plt.figure("Shannon's limit", figsize=(8,6))
plt.xlim(-2-5*dx,11+5*dx)
plt.ylim(-10*dy,7)
plt.axis('off')

ruler(X(-2),X(10),right_inc=15)
ruler(Y(0),Y(6),right_inc=15)


# Ticking utilities and ticks on axes
def vtick(x,l=0.01):
    return seg((x,-l),(x,l))
def htick(y,l=0.01):
    return seg((-l,y),(l,y))
#
for a in [-2,2,4,6,8,10]:
    vtick(a,0.6*dx)
for b in [1,2,3,4,5,6]:
    htick(b,dy)
#
vtick(-1.59,0.6*dx)

# Dashed segments
seg((R[0],-8*dy),R,dashing='--')
seg((B[0],-8*dy),B,dashing='--')
seg((L[0],-8*dy),(L[0],0),dashing='--')

seg(shannon_point(1),B,dashing='--',color='blue')
seg(shannon_point(1/3),R,dashing='--',color='red')
seg((shannon_limit(1/3),-2*dy),shannon_point(1/3),dashing='--')

# Drawing of Shannon's limiting curve
plot(shannon_limit(s),s,color='blue')

# Bullets
bullet((shannon_limit(1/3),1/3),color='red')
bullet(shannon_point(1),color='blue')
bullet((6.26,1/3),color='red'); bullet((9.6,1),color='blue')

# Lables
lable((9.6,1),r'$B$', dx=dx,dy=dy,fs=16, color='blue')
lable((6.26,1/3),r'$R$', dx=dx,dy=dy,fs=16, color='red')

lable(X(-2),r'$-2$',dx=-6*dx,dy=-5*dy,fs=16)
lable(X(2),r'$2$',dx=-1.5*dx,dy=-5*dy,fs=16)
lable(X(4),r'$4$',dx=-1.5*dx,dy=-5*dy,fs=16)
lable(X(6),r'$6$',dx=-1.5*dx,dy=-5*dy,fs=16)
lable(X(8),r'$8$',dx=-1.5*dx,dy=-5*dy,fs=16)
lable(X(10),r'$10$',dx=-2*dx,dy=-5*dy,fs=16)

lable((shannon_limit(1/3),-4*dy),r'$-1.08$',dx=-2*dx,dy=0,fs=12)
lable((L[0],-8*dy),r'$-1.59$',dx=dx,dy=0,fs=12)
lable((R[0],-8*dy),r'$6.26$',dx=dx,dy=0,fs=12)
lable((B[0],-8*dy),r'$9.59$',dx=-9.5*dx,dy=0,fs=12)


for j in range(1,7):
    lable(Y(j),str(j),dx=-6*dx,dy=0,fs=16)

lable((11.5,0),r'$\rho$',dx=0,dy=-3*dy,fs=16)
lable((0,6.6),r'$\eta$',dx=2*dx,dy=0,fs=16)

    
plt.show()